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Abstract 

R versions of the array manipulation functions of APL are presented. We 
do not translate the system functions or other parts of the runtime. Also, the 
current version has does not have the nested arrays of APL-2. 
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the R, C, and C++ source code. 

1 Introduction 

APL was introduced by Iverson (1962). It is an array language, with many functions 
to manipulate multidimensional arrays. R also has multidimensional arrays, but 
not as many functions to work with them. 

In R there are no scalars, there are vectors of length one. For a vector x in R we 
have dim(x) equal to NULL and length(x) > 0. For an array, including a matrix, we 
have length (dim (x)) > 0. APL is an array language, which means everything is an 
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array. For each array both the shape pA and the rank ppA are defined. Scalars are 
arrays with shape equal to one, vectors are arrays with rank equal to one. 

If you want to evaluate APL expressions using a traditional APL virtual keyboard, 
we recommend the nice webpage at ngn.github.io/apl/web/index.html. EliStudio 
at fastarray.appspot.com/default.html is essentially an APL interpreter running in 
a Qt GUI, using ascii symbols and symbol-pairs to replace traditional APL symbols 
(Chen and Ching (2013)). Eli does not have nested arrays. It does have ecc, which 
compiles eli to C. 

In 1994 one of us coded most APL array operations in XLISP-STAT. The code is still 
available at gifi.stat.ucla.edu/apl. 

There are some important differences between the R and Lisp versions, because 
Lisp and APL both have C’s row-major ordering, while R (like Matlab) has Fortran’s 
column-major ordering in the array layout. Our R version of APL uses column- 
major ordering. By slightly changing the two basic building blocks of our code, 
the aplDecodeO and aplEncodeQ functions, it would be easy to choose between 
row-major and column-major layouts. But this would make it more complicated to 
use the code with the rest of R. 

Because of layout, the two arrays 3 3 3pr27 and array (1: 27, rep (3,3)) are different. 
But what is really helpful in linking the two environments is that ,3 3 3pr27 and 
as .vector(array(1:27,rep(3,3)), which both ravel the array to a vector, give the 
same result, the vector a.27 or 1:27. This is, of course, because ravelling an array 
is the inverse of reshaping a vector. 

Most of the functions in R are written with arrays of numbers in mind. Most of 
them will work for array with elements of type logical, and quite a few of them 
will also work for arrays of type character. We have to keep in mind, however, 
that APL and R treat character arrays quite differently. In R we have length("aa") 
equal to 1, because "aa" is a vector with as its single element the string "aa". R 
has no primitive character type, characters are just strings which happen to have 
only one character in them. In APL strings themselves are vectors of characters, 
and paa is 2. In R we can say a<-array ("aa", c(2,2,2)), but in APL this gives a 
domain error. In APL we can say 2 2 2p“aa”, which gives the same result as 2 2 
2p“a” or 2 2 2p‘a\ 

In this version of the code we have not implemented the nested arrays of APL-2. 
Nesting gives every array A not just a shape pA and a rank ppA, but also a depth. 
The depth of an array of numbers or characters is one, the depth of a nested array 
is the maximum depth of its elements. 

There are many dialects of APL, and quite a few languages derived from APL, such 
as A+ and J. As a standard for APL-I we use Helzer (1989). 
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2 Core 


The basic engine behind the APL array manipulation functions is the pair 
aplDecodeQ and aplEncodeO. They make it easy to go back and forth be¬ 
tween multidimensional arrays and the one-dimensional vectors that store their 
elements. 

Suppose a is the array 


## , , 1 
## 


## 


[,1] 

[, 2] 

[, 3] 

## 

[1J 

1 

3 

5 

## 

[2,] 

2 

4 

6 

## 





## 

r 

) 



## 





## 


[,1] 

[, 2] 

[, 3] 

## 

[1J 

7 

9 

11 

## 

[2,] 

8 

10 

12 

## 





## 

} } ^ 

5 



## 





## 


[,1] 

[, 2] 

[, 3] 

## 

[1J 

13 

15 

17 

## 

[2,] 

14 

16 

18 

## 





## 

J > 4 

1 



## 





## 


[,1] 

[, 2] 

[, 3] 

## 

[1J 

19 

21 

23 

## 

[2,] 

20 

22 

24 


with shape equal to 
## [1] 2 3 4 
and rank 
## [1] 3 

Of course the length of a is the product of the elements of its shape. If r is an 
integer between 1 and length (a) then aplEncode (r, aplShape (a)) returns the index 
vector k such that element with indices k of a is a[r] . 
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aplSelect(a, aplEncode (1, aplShape (a)), drop = TRUE) 


## [1] 1 

aplSelect(a, aplEncode (10, aplShape(a) ), drop = TRUE) 

## [1] 10 

aplSelect(a, aplEncode (24, aplShape (a)), drop = TRUE) 

## [1] 24 

If k is an admissible index vector then aplDecode (k, aplShape (a)) returns an integer 
such that element with indices k of a is a [aplDecode (k, dims) ] 

a[aplDecode (c(l,l,l), aplShape(a) )] 

## [ 1 ] 1 

a[aplDecode (c(2,2,2), aplShape (a))] 

## [ 1 ] 10 

a[aplDecode (c(2,3,4), aplShape(a) )] 

## [1] 24 

Note that aplDecodeO and aplEncodeO are inverses of each other because 
r <- 2 

r == aplDecode(aplEncode (r, aplShape(a) ), aplShape(a)) 

## [1] TRUE 
k <- c ( 1 , 2 ,3) 

all(k == aplEncode(aplDecode (k, aplShape (a)), aplShape (a))) 

## [1] TRUE 

and this is true for all admissible r and k. 
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3 Functions 


3.1 Compress 

Compress (IBM (1988) , p. 91—92) is defined as the Replicate operator L/R in the 
special case that L is binary. So look under Replicate for the definition. 


3.2 Decode 

The decode dyadic operator L_LR is also known as base value (Helzer (1989), p. 17-21). 
If L is scalar and R is a vector, then L_LR is the polynomial rqr ” 1 ^ 1 +r 2 x m ~ 2 + • • • + r rn 
evaluated at L. This means that if the r* are nonnegative integers less than L, then 
L_LR gives the base-10 equivalent of the base-L number R. 

Normally, however, and in our R implementation, the arguments L and R are 
vectors of the same length. This is also because for scalar L the expression L_LR 
is expanded to ((pR)pL)_LR If L and R are vectors of the same length then decode 
returns the index of element A [R] in an array A with dim(A)=L. 

Obviously the R implementation, which uses colum-major ordering, will give re¬ 
sults different from the APL implementation. In APL the expression 3 3 3_L1 2 3 
evaluates to 18, while aplDecode(l: 3, rep(3,3)) gives 22. Also note the order of 
the arguments is interchanged. Also note that if x and y are of length m, and y 
has all elements equal to z, then aplDecode(x, reply, length(x))) is the value of 
the polynomial 

p(u) — 1 + (xi — l)w° + (x 2 - l)^ 1 4-b (x m - l)u m L 

evaluated at u — z. 

for (k in 1:3) for (j in 1:3) for (i in 1:3) 
print (aplDecode (c(i,j,k), c(3,3,3))) 


## 

[1] 

1 

## 

[1] 

2 

## 

[1] 

3 

## 

[1] 

4 

## 

[1] 

5 

## 

[1] 

6 

## 

[1] 

7 

## 

[1] 

8 

## 

[1] 

9 

## 

[1] 

10 

## 

[1] 

11 

## 

[1] 

12 
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## [1] 13 
## [1] 14 
## [1] 15 
## [1] 16 
## [1] 17 
## [1] 18 
## [1] 19 
## [ 1 ] 20 
## [ 1 ] 21 
## [ 1 ] 22 
## [1] 23 
## [1] 24 
## [1] 25 
## [1] 26 
## [1] 27 


3.3 Drop 

See Helzer (1989), p. 49—51. If L is a positive integer and R is a vector then L'lR 
drops the first L elements. If L is a negative integer the last L elements are 
dropped. If R is an array, then L must have ppR elements. Depending on whether 
the elements or L are positive or negative, the L first or last slices of the dimension 
are dropped. 

Our R implementation again interchanges the two arguments. Note that in R 
the default on subsetting arrays is drop = TRUE. In our implementations the de¬ 
fault is always drop = FALSE. Note that more general subsetting can be done with 
aplSelect(). 

So for a vector 
aplDrop(l : 10,3) 


## [1] 4 5 6 7 8 9 10 


aplDrop(l : 10,-3) 


## [1] 1 2 3 4 5 6 7 


and for an array 

aplDrop(a, c(l , 0,1) , drop = TRUE) 
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## [,1] C,2] [, 3] 

## [1,] 8 14 20 

## [2,] 10 16 22 

## [3,] 12 18 24 

aplDrop(a, c(-l,-l,0) , drop = TRUE) 

## [,1] [,2] [,3] [,4] 

## [1J 1 7 13 19 

## [2,] 3 9 15 21 

3.4 Encode 

Encode, also known as representation is the inverse of decode (Helzer (1989), p. 17— 
21). LTR takes a radix vector L and a number R and returns the array indices 
corresponding to cell L in an array with shape pA. 

for (i in 1:27) 

print (aplEncode (i, c(3,3,3))) 


## [ 1 ] 1 1 1 

## [ 1 ] 2 1 1 

## [1] 3 1 1 

## [ 1 ] 1 2 1 

## [ 1 ] 2 2 1 

## [1] 3 2 1 

## [1] 13 1 
## [1] 2 3 1 

## [1] 3 3 1 

## [1] 1 1 2 

## [1] 2 1 2 

## [1] 3 1 2 

## [1] 12 2 
## [ 1 ] 2 2 2 

## [1] 3 2 2 

## [ 1 ] 13 2 
## [1] 2 3 2 

## [1] 3 3 2 

## [1] 1 1 3 

## [1] 2 1 3 

## [1] 3 1 3 

## [1] 12 3 
## [1] 2 2 3 
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## [1] 3 2 3 

## [1] 13 3 

## [1] 2 3 3 

## [1] 3 3 3 

3.5 Expand 

Expand (Helzer (1989), p. 64—66) L\R replaces slices of a vector or an array by 
zeroes. In APL we use a boolean vector for L, and an array for R. If R is a vector 
then the number of elements of L equal to one (i.e. TRUE) should be equal to the 
length of R. Then L\R produce a vector of length pL with the elements of R in 
the places where L is one, and zeroes elsewhere. In our function aplExpandO we 
again reverse the order of the arguments. The second argument in the R verson 
can be both logical or binary. 

aplExpand(l : 3 , c(l,0,0,0,1,1)) 

## [1] 1 0 0 0 2 3 

For an array expand takes an optional axis argument. For a matrix, for example, 
axis=l expands rows (i.e. inserts row of zeroes at specfied places), while axis=2 
expands columns. This generalizes to multidimensional arrays, where there are 
just more axis to consider, but any one of them can be expanded. 

aplExpand(matrix (1,2,3) , c(l,0,0,l), axis = 1) 

## [, 1] [, 2] [, 3] 

## [ 1 ,] 1 1 1 

## [2,] 000 

## [3,] 0 0 0 

## [4,] 1 1 1 

aplExpand(matrix(l ,2,3) , c(l,l,0,1,0) , axis = 2) 

## [,1] [, 2] [, 3] [, 4] [, 5] 

## [ 1 ,] 11010 

## [ 2 ,] 11010 

3.6 Get 

There is no APL function corresponding to get. We wrote aplGetO as a simple 
utility to get an element out of an array, 
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a <- array(l:24, c(2,3,4)) 
aplGet (a, c(2,2,2)) 

## [ 1 ] 10 

aplGet (a, aplEncode (14, c(2,3,4))) 

## [1] 14 

3.7 Inner Product 

APL has a generalized inner product (Helzer (1989), p. 100-103), which we write as 
Lf.gR, where f and g are any scalar dyadic functions. In obvious notation the f,g 
inner product of two vector x and y of length n is f{g(xi,yi),g(x 2 ,y 2 ),''' , g{x n , y n ))- 
If L and R are arrays, then the last item of pL must be equal to the first item of 
pR and the generalized inner product operates along that common dimension. The 
default for g is multiplicaton, and for f is addition, leading to the ordinary inner 
product and to matrix multiplication. 

Of course there are many examples we can give. 


x <- matrix (1:12,4,3) 
y <- matrix (1:12,3,4) 
aplInnerProduct (x, y) 


## 


[,1] 

[, 2] 

[, 3] 

[ > 4] 

## 

[1J 

38 

83 

128 

173 

## 

[2,] 

44 

98 

152 

206 

## 

[3,] 

50 

113 

176 

239 

## 

[4,] 

56 

128 

200 

272 

h < 

- function 

(x, 

y) ifelse 

aplInnerProduct ( 

x, y, 

h, " 

## 


[,1] 

[, 2] 

[, 3] 

[ > 4] 

## 

[1J 

1 

1 

1 

0 

## 

[2,] 

0 

0 

0 

0 

## 

[3,] 

0 

0 

0 

0 

## 

[4,] 

0 

1 

1 

1 


y, i, o) 
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a <- array (1:24, c(2,3,4)) 
b <- rep (1, 4) 
aplInnerProduct (a, b) 

## [,1] [, 2] [,3] 

## [1,] 40 48 56 

## [2,] 44 52 60 


3.8 Join 

Join catenates two arrays of the same rank, creating another larger array of that 
rank (Helzer (1989), p. 104—110). If L and R are vectors L,R simply concatenates, 
same as c() in R. If L and R are arrays they can be catenated along axis k if (pL)[-k] 
= (pR)[-k]. Thus matrices with the same number of columns can be stacked on top 
of each other, and matrices with the same number of rows can be stacked next to 
each other. Same for higher-dimensional arrays. 

In APL there are some special rules for handling scalar arguments (they get re¬ 
shaped accordingly first), and even for fractional axis parameters, which can be 
used for lamination of arrays with different rank. We have not implemented 
these more complicated optiosn in our function aplJoinO. 

x<-l : 3 
y<-3:l 
aplJoin(x,y) 

## [1] 1 2 3 3 2 1 

x <- matrix (1:12, 3, 4) 
y <- matrix (1:8, 2, 4) 
aplJoin (x, y, axis = 1) 


## [, 1] [, 2] [, 3] 

## [1,] 1 4 7 

## [ 2 ,] 258 

## [3,] 3 6 9 

## [4,] 1 3 5 

## [5,] 246 


[,4] 

10 

11 

12 

7 

8 


a <- array (1:24, c(2, 3, 4)) 
b <- array (1:30, c(2, 3, 5)) 
dim(aplJoin(a, b, axis = 3)) 


## [1] 2 3 9 


11 


3.9 Member Of 


The member of function in APL LER returns a binary array with the shape of L. 
Array R can be of any shape. If an element of L occurs in R then the corresponding 
element in LER is one, otherwise it is zero. Our function aplMemberOf () returns 
an array or vector of the same dimension or length as the first argument. The 
test for equality is a simple ==, so integers and doubles can be compared. 

a <- array (1:24, c(2,3,4)) 
aplMemberOf (a, c(l , 2 , 15 ,25) ) 


## , , 1 

## 

## [, 1] [, 2] [, 3] 

## [1,] 1 0 0 

## [2,] 1 0 0 

## 

## , , 2 
## 

## [, 1] [, 2] [, 3] 

## [1,] 000 

## [2,] 000 

## 

## , , 3 

## 

## [,1] [, 2] [, 3] 

## [1,] 0 1 0 

## [2,] 000 

## 

## , , 4 
## 

## [, 1] [,2] [, 3] 

## [1,] 000 

## [2,] 000 


3.10 Outer Product 

The APL outer product L °.f R is the same as the outer product in R (Helzer (1989), 
p. 147—149). Both L and R can be arbitrary arrays, and the result is an array with 
shape (pL),pR, formed by applying the dyadic function f to each combination of 
elements of L and R. Thus our aplOuterProduct() just calls the R function outer() 
on its arguments. 
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x <- matrix (1:4, 2,2) 
y <- matrix (1:4, 2,2) 
aplOuterProduct (x, y, "+") 


## , , 1, 1 

## 

## [, 1 ] [, 2 ] 
## [1,] 2 4 

## [2,] 3 5 

## 

## , , 2 , 1 

## 

## [, 1 ] [, 2 ] 

## [1,] 3 5 

## [2,] 4 6 

## 

## , , 1 , 2 
## 

## [, 1 ] [, 2 ] 

## [1,] 4 6 

## [2,] 5 7 

## 

## , , 2 , 2 

## 

## C, 1] [, 2] 

## [1,] 5 7 

## [2,] 6 8 


3.11 Ravel 

The APL function ravel reshapes its argument into a vector. So ,R is identical to 
as.vectorO in R, and aplRavelO is just a call to as.vector(). 

aplRavel(array (1 : 24, c(2,3,4))) 

## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 

## [24] 24 


3.12 Rank 

There is no APL function rank, but we just implemented aplRankO as a convenient 
shorthand for ppR. In R we just call length (dim ()). 
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aplRank(array (1: 24, c(2,3, 4))) 


## [1] 3 

3.13 Reduce 

It is time for a quote from Helzer (1989), p. 165. “Summing a list of numbers is a 
common programming task. In APL this, and much more, is accomplished using 
the reduction operator. In APL terminology an operator modifies the action of 
a function. Standard APL has four operators, inner product (f.g), outer product 
(°.f), reduction (f/), and scan (f\). Inner product creates a new dyadic function 
out of two scalar dyadic functions f and g. Outer product creates a new dyadic 
functionout of a scalar dyadic function f. Reduction and scan both create a new 
monadic function out of a scalar dyadic function f.” In APL there are two operators 
which reduces along the first axis and /, which reduces along the last axis. 
Reducing along the k-th axis is f/[k]. 

Our R version aplReduceO of reduce takes three arguments: the array, the axes 
to reduce over, and the dyadic function used for reduction (by default addition). 
Note that the second argument can be a vector with more than one axis, in which 
case we reduce over all those axes. 

a <- array(l:24, c(2,3,4)) 
aplReduce (a, c(l,2), max) 

## [1] 6 12 18 24 

aplReduce (a, 1, function (x, y) ifelse (x > y, x, y)) 


## 


[,1] 

[, 2] 

[, 3] 

[>4] 

## 

[1J 

2 

8 

14 

20 

## 

[2,] 

4 

10 

16 

22 

## 

[3,] 

6 

12 

18 

24 


3.14 Replicate 

If L and R in replicate L/R are vectors of the same length, then the result is a 
vector of length +/L in which element r[i] is repeated l[i\ times. If L is binary 
then some elements will be deleted, and replicate is called compress. 
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aplReplicate (1 :3,c(3,1,3)) 


## [ 1 ] 1112333 
aplReplicate (1: 10 , rep (c(0,1) ,5) ) 


## [1] 2 4 6 8 10 

This can be extended to higher dimensional arrays in the usual way by specifying 
the axis along which to replicate. By default, as in APL, we select the last axis. 
In APL L/R is used to replicate along the first axis. Our function aplReplicate() 
covers both cases. 

a <- array (1:24,c(2,3,4)) 

aplReplicate(aplReplicate (a, c(2,2), 1), c(0,2,0), 2) 


## 

, , 1 



## 




## 


[,1] 

[, 2] 

## 

[1J 

3 

3 

## 

[2,] 

3 

3 

## 

[3,] 

4 

4 

## 

[4,] 

4 

4 

## 




## 

, , 2 



## 




## 


[,1] 

[, 2] 

## 

[1J 

9 

9 

## 

[2,] 

9 

9 

## 

[3,] 

10 

10 

## 

[4,] 

10 

10 

## 




## 

, , 3 



## 




## 


[,1] 

[, 2] 

## 

[1J 

15 

15 

## 

[2,] 

15 

15 

## 

[3,] 

16 

16 

## 

[4,] 

16 

16 

## 




## 

, , 4 



## 
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## [, 1 ] [, 2 ] 

## [ 1 ,] 21 21 

## [ 2 ,] 21 21 

## [3,] 22 22 

## [4,] 22 22 

aplReshape (c(l,2), c(2,2,2)) 


## , , 1 

## 

## [, 1 ] [, 2 ] 

## [ 1 ,] 1 1 

## [ 2 ,] 2 2 

## 

##,, 2 
## 

## [, 1 ] [. 2 ] 

## [ 1 ,] 1 1 

## [ 2 ,] 2 2 


aplReshape (array(l:24, c(2,3,4)), c(2,2)) 


## [, 1 ] [, 2 ] 

## [1,] 1 3 

## [2,] 2 4 


3.15 Rotate 

Rotate (Helzer (1989), p. 191 — 193) cyclically shifts the elements of a vector or array 
dimension. In APL we write either L(J)R or L©R, depending on whether rotation 
is along the first or the last axis. 

For vectors aplRotate(x, k) with positive k takes the first k elements of x and 
moves them to the end, with negative k it takes the last k elements and moves 
them to the front. 

aplRotate(l : 6, 2) 

## [1] 3 4 5 6 1 2 
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aplRotate(l : 6, -2) 


## [1] 5 6 1 2 3 4 


For a matrix we can shift rows or columns. For a matrix R the expression L(J)R, 
or aplRotate (R, L, axis = 2), shifts the elements of rows, with within each row 
possibly a different shift. Thus the number of elements of L must be the number 
of rows of R. If L is a scalar, then it is basically first blown up into a vector, and 
the same shift L is applied to each row (which means that a column is shifted). If 
axis = 1 then elements within columns are shifted. By default the argument axis 
equals aplRank(a). 

a <- cbind(l:5, matrix (0,5,4)) 
aplRotate (a, 2) 


## [,1] [, 2] [, 3] 

## [ 1 ,] 000 

## [ 2 ,] 000 

## [3,] 0 0 0 

## [4,] 000 

## [5,] 000 


[,4] [, 5] 
1 0 

2 0 

3 0 

4 0 

5 0 


aplRotate (a, -c(0, 1 ,2,3,4)) 


## 

## [ 1 ,] 
## [ 2 ,] 
## [3,] 
## [4,] 
## [5,] 


[,1] [, 2] [,3] [, 4] [, 5] 
1 0 0 0 0 

0 2 0 0 0 

0 0 3 0 0 

0 0 0 4 0 

0 0 0 0 5 


aplRotate (a, 2, 1) 


## 

## [ 1 ,] 
## [ 2 ,] 
## [3,] 
## [4,] 
## [5,] 


[, 1 ] [, 2 ] 

3 0 

4 0 

5 0 

1 0 

2 0 


[, 4] [, 5] 
0 0 

0 0 

0 0 

0 0 

0 0 


[,3] 

0 

0 

0 

0 

0 
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For arrays matters become more complicated. To rotate R along k we need dim(L) 
to be dim(R) [—k]. Then each element in L defines a slice of R that is a vector of 
length dim(R)[k]. That vector then gets rotated using the positive or negative 
value of the corresponding element of L. 

a <- array(l:24, c(2,3,4)) 
aplRotate (a, 1, 1) 


## , , 1 
## 

## C, 1] [, 2] [,3] 

## [1,] 2 4 6 

## [2,] 1 3 5 

## 

## , , 2 
## 

## [, 1] [, 2] [,3] 

## [ 1 ,] 8 10 12 

## [2,] 7 9 11 

## 

##,, 3 
## 

## [,1] [,2] [,3] 

## [1,] 14 16 18 

## [2,] 13 15 17 

## 

## , , 4 
## 

## C, 1] [, 2] [, 3] 

## [1,] 20 22 24 

## [2,] 19 21 23 

b <- matrix (c(0,0,0,l,l 
aplRotate (a, b, 3) 


1) , 2 , 3, byrow = TRUE) 


## 

## 

## 

, , 1 

[,1] 

[,2] 

[, 3] 

## 

[1J 

1 

3 

5 

## 

[2,] 

8 

10 

12 

## 

## 

## 

, , 2 
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## 


[,1] 

[,2] 

[, 3] 

## 

[1J 

7 

9 

11 

## 

[2,] 

14 

16 

18 

## 





## 

, , 3 




## 





## 


[,1] 

[, 2] 

[, 3] 

## 

[1J 

13 

15 

17 

## 

[2,] 

20 

22 

24 

## 





## 

, , 4 




## 





## 


[,1] 

[, 2] 

[, 3] 

## 

[1J 

19 

21 

23 

## 

[2,] 

2 

4 

6 

3.16 Scan 




The APL operator scan does not change dimension (Helzer (1989), p. 195—197). The 
expression f\R for a vector R can be defined in terms of reduction. It produces 
the vector with elements f/R[l] , f/R[l 2] , f/R[l 2 3], .... Thus it generalizes R 
funcions such as cumsumO, cumprodO, and so on. For a matrix f\R scans the rows 
and f\[l]R or f\R scans the columns. For arrays we again need an axis to expand 
along. 

a<-matrix(l :9,3,3) 
aplScan(a, 1, "+") 


## [,1] [, 2] [, 3] 

## [1J 1 4 7 

## [2,] 3 9 15 

## [3,] 6 15 24 

aplScan(a, f = min) 

## [,1] [. 2] [,3] 

## [1J 1 1 1 

## [ 2 ,] 222 

## [3,] 3 3 3 


a<-array(l:24,c(2,3,4)) 
aplScan(a, 3, "*") 
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## 

> > 

1 



## 





## 


[,1] 

[,2] 

[, 3] 

## 

[1J 

1 

3 

5 

## 

[2,] 

2 

4 

6 

## 





## 

J J 

2 



## 





## 


[,1] 

[, 2] 

[,3] 

## 

[1J 

7 

27 

55 

## 

[2,] 

16 

40 

72 

## 





## 

> > 

3 



## 





## 


[,1] 

[, 2] 

[, 3] 

## 

[1J 

91 

405 

935 

## 

[2,] 

224 

640 

1296 

## 





## 

> > 

4 



## 





## 


[,1] 

[,2] 

[,3] 

## 

[1J 

1729 

8505 

21505 

## 

[2,] 

4480 

14080 

31104 


3.17 Select 

Select is not an APL function. Our function aplSelectO has an array a as its first 
argument, and a list of aplRank(a) vectors of integers as it second element. It 
then selects the corresponding rows, columns, slices, and so on. 

a <- array(l:24, c(2,3,4)) 
aplSelect (a, list(l,c(l,2),c(3,4))) 

## C, 1] [, 2] 

## [1,] 13 19 
## [2,] 15 21 

aplSelect(a, list (1 , c(l ,2) , c(3 ,4) ), drop = FALSE) 


## , , 1 
## 

## C, 1] [, 2] 

## [1,] 13 15 
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## 

##,, 2 
## 

## C, 1] [, 2] 

## [1,] 19 21 


3.18 Set 

There is no APL function corresponding with set. We wrote aplSetO as a simple 
utility to set an element of an array to a value. 

a <- array(l:12, c(2,3,2)) 
aplSet (a, 11, c(2,2,2)) 


## , , 1 
## 

## C, 1] [, 2] [, 3] 

## [1,] 1 3 5 

## [ 2 ,] 246 

## 

## , , 2 
## 

## C, 1] [, 2] [, 3] 

## [1,] 7 9 11 

## [ 2 ,] 8 11 12 


3.19 Shape 

Shape is the monadic version of p, while reshape is the dyadic version. Shape gives 
the dimensions of an array, an reshape modifies them. Of course aplShapeO is 
just a wrapper for the basic R function dim(). 

a <- array(l:12, c(2,3,2)) 
aplShape(a) 


## [1] 2 3 2 
aplShape(aplShape (a)) 


## [1] 3 
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3.20 Take 


Take LtR extracts items from vectors and arrays, in the same way as drop iTR 
deletes items. If R is a vector then L items from the beginning of R are taken if 
L is positive, and L items from the end if L is negative. If R is an array, then L 
must have ppR elements. Depending on whether the elements or L are positive or 
negative, the L first or last parts of the dimension are taken 

Our R implementation again interchanges the two arguments. Note that in base 
R the default on subsetting arrays is drop = TRUE. In our implementations the 
default is always drop = FALSE. Both aplTakeO and aplDropO are implemented 
by constructing a suitable list of index vectors for aplSelect(). Note that more 
general subsetting can be done with aplSelect(). Note that 

So for a vector 
aplTake(l : 10,3) 


## [1] 12 3 
aplTakeCl: 10,-3) 


## [1] 8 9 10 

and for an array 

a <- array(l:24, c(2,3,4)) 
aplTake(a, c(2,3,2), drop = TRUE) 


## 

> 

, 1 




## 






## 



[,1] 

[,2] 

[, 3] 

## 

[1 

J 

1 

3 

5 

## 

[2 

J 

2 

4 

6 

## 






## 

> 

, 2 




## 






## 



[,1] 

[, 2] 

[, 3] 

## 

[1 

J 

7 

9 

11 

## 

[2 

J 

8 

10 

12 


22 


aplTake(a, c(2,-2,l), drop = TRUE) 


## [, 1 ] [, 2 ] 

## [1,] 3 5 

## [2,] 4 6 

3.21 Transpose 

APL has both a monadic (S)R and a dyadic L(S)R transpose. This APL transpose has a 
somewhat tortuous relationship with apermO in R. 

The monadic aplTranspose(a) and aperm(a) are always the same, they reverse the 
order of the dimensions. 

If x is a permutation of 1:aplRank(a), then aperm(a,x) is actually equal to 
aplTranspose (a, order (x)) For permutations we could consequently define 
aplTranspose(a,x) simply as aperm(a,order(x)) (which would undoubtedly be 
more efficient as well). 

a <- array(l:24, c(2, 3, 4)) 
aplTranspose (a) 


## , , 1 
## 


## 


[,1] 

[, 2] 

[, 3] 

## 

[1J 

1 

3 

5 

## 

[2,] 

7 

9 

11 

## 

[3,] 

13 

15 

17 

## 

[4,] 

19 

21 

23 

## 





## 

, , 2 




## 





## 


[,1] 

[, 2] 

[, 3] 

## 

[1J 

2 

4 

6 

## 

[2,] 

8 

10 

12 

## 

[3,] 

14 

16 

18 

## 

[4,] 

20 

22 

24 


aplTranspose (a, c(2,l,3)) 


## , , 1 
## 

## [, 1 ] [, 2 ] 
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## [ 1 ,] 1 2 

## [2,] 3 4 

## [3,] 5 6 

## 

##,, 2 
## 

## C, 1] [, 2] 

## [1,] 7 8 

## [2,] 9 10 

## [3,] 11 12 

## 

##,, 3 
## 

## [, 1 ] [, 2 ] 

## [1,] 13 14 

## [2,] 15 16 

## [3,] 17 18 

## 

## , , 4 
## 

## [, 1 ] [, 2 ] 

## [1,] 19 20 

## [ 2 ,] 21 22 

## [3,] 23 24 

If x is not a permutation, then aperm(a,x) is undefined, but aplTranspose(a,x) 
can still be defined in some cases. If x has aplRank(a) elements equal to one of l:m, 
with each of l:m occurring a least once, then aplTranspose(a,x) has rank m. If an 
integer between 1 and m occurs more than once, then the corresponding diagonal 
of a is selected. 

a 


## 

, , 1 



## 




## 

[,1] 

[, 2] 

[, 3] 

## 

[1J 1 

3 

5 

## 

[2,] 2 

4 

6 

## 




## 

, , 2 



## 




## 

[,1] 

[, 2] 

[, 3] 

## 

[1J 7 

9 

11 

## 

[2,] 8 

10 

12 


24 



## 





## 

, , 3 




## 





## 


[,i] 

[, 2] 

[, 3] 

## 

[1J 

13 

15 

17 

## 

[2,] 

14 

16 

18 

## 





## 

, , 4 




## 





## 


[,1] 

[, 2] 

[, 3] 

## 

[1J 

19 

21 

23 

## 

[2,] 

20 

22 

24 

aplTranspose 

(a, 

c (2 , 2 , 1) ) 

## 


[,1] 

[, 2] 


## 

[1J 

1 

4 


## 

[2,] 

7 

10 


## 

[3,] 

13 

16 


## 

[4,] 

19 

22 


4 

Implementations 


4.1 Six Versions 

We have implemented the APL functions in six different ways 

1. As pure R. 

2. As R, with only decode and encode in C using .C(). 

3. As C, using the ,C() interface for decode, encode, transpose, select, reduce, 
scan, and inner product. 

4. As C, using the .Call() interface for decode, encode, transpose, select, reduce, 
scan, and inner product. 

5. As C, using the .Call() interface, for transpose, select, reduce, scan, and 
inner product, with decode and encode inlined using .C(). 

6. As C using Repp for decode, encode, transpose, select, reduce, scan, and inner 
product. 

It must be emphasized that the Repp interface was written five years ago, with a 
very early version of Repp, that uses only a tiny subset of the possibilities offered 
by newer versions. We are sure a great deal of improvement is possible there. 
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Also for reduce, scan, and inner product some of our arguments are functions. 
In the ,C() implementation we use the old Call_R() interface, dating back to the 
Blue Book (Becker, Chambers, and Wilks (1988)), to handle function pointers. 

Note that inlining decode and encode makes sense, since they are called so many 
times in all functions, but the inline keyword is handled only as a hint to the 
compiler. It does not guarantee actual inlining of the function. 

The computations in the body of the paper use the Call() interface. The different 
implementations calling C routines use different shared libraries and different 
glue routines in R. 


4.2 Timing 

We compare running time of the six implementations. The parameters we use are 
repNum <- 10 

timeArr <- array (NA, c (repNum, 3,3, 6)) 
a<-array(l : 10000, c(10, 10, 100)) 
b<-array(l : 10000, c( 100, 10,10)) 
c<-array(l : 100000, rep(10, 5)) 
d<-array(l : 100000, rep(10, 5)) 
x<-list (1: 5 ,1: 5 , 1 : 5 , 1 : 5 ,1: 5) 


and the results are 


## , , InnerProduct 
## 

## R R+.C .C .Call Inline Repp 
## user 9.7987 10.8804 0.8723 0.3620 0.2348 0.7785 
## system 0.0372 0.0635 0.0491 0.0046 0.0030 0.0100 
## elapsed 9.9047 11.1021 0.9257 0.3687 0.2398 0.8048 


## 




## 

, , Reduce 


## 




## 


R 

R+.C 

## 

user 

1.6892 

1.0227 

## 

system 

0.0086 

0.0071 

## 

elapsed 

1.7131 

1.0435 

## 




## 

, , Select 


## 




## 


R 

R+.C 

## 

user 

0.0624 

0.0497 

## 

system 

0.0002 

0.0005 

## 

elapsed 0.0627 

0.0563 


.C .Call Inline Repp 
0.0464 0.0476 0.0322 0.0732 
0.0034 0.0016 0.0010 0.0025 
0.0498 0.0495 0.0335 0.0767 


.C .Call 
8e-04 8e-04 
0e+00 0e+00 
7e-04 9e-04 


Inline Repp 
5e-04 0.0024 
0e+00 0.0001 
6e-04 0.0024 
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5 Appendix: Code 

5.1 R Code 


# select 

aplSelect <- function(a, x, drop = FALSE) { 
sa <- aplShape(a) 
ra <- aplRank(a) 
sz <- sapply(x, length) 
z <- array(0, sz) 
nz <- prod(sz) 
for (i in l:nz) { 

ivec <- aplEncode(i, sz) 
jvec <- vectorO 
for (j in lira) 

jvec <- c(jvec, x [ [j] ] [ivec [j] ]) 
z[i] <- a[aplDecode(jvec, sa)] 

> 

if (drop) 

return(drop(z) ) 
else 

return(z) 


# drop 

aplDrop <- function(a, x, drop = FALSE) { 
sa <- aplShape(a) 
ra <- aplRank(a) 
y <- as.list(rep(0 , ra)) 
for (i in lira) { 
ss <- sa[i] 
xx <- x[i] 
sx <- ss + xx 
if (xx >= 0) 

y[[i]] <- (xx + 1) : ss 
if (xx < 0) 
y C [i]] <- 1 :sx 

> 

return(aplSelect (a, y, drop)) 
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# take 


aplTake <- function(a, x, drop = FALSE) { 
sa <- aplShape(a) 
ra <- aplRank(a) 
y <- as.list(rep(0 , ra)) 
for (i in l:ra) { 
ss <- sa[i] 
xx <- x[i] 
sx <- ss + xx 
if (xx > 0) 

y[[i]] <- l:xx 
if (xx < 0) 

y[[i]] <- (sx + 1) : ss 

} 

return(aplSelect (a, y, drop)) 


# reduce vector 

aplRDV <- function(x, f = "+") { 
if (length(x) == 0) 
return(x) 
s <- x[l] 

if (length(x) == 1) 
return(s) 

for (i in 2: length (x)) 

s <- match.fun (f)(s, x[i]) 
return(s) 


# scan vector 

aplSCV <- function(x, f = "+") { 
if (length(x) <= 1) 
return(x) 

return(sapply(l:length(x) , function(i) 
aplRDV (x[l :i], f))) 

> 

# inner product vector 

aplIPV <- function(x, y, f = g = "+") { 

if (length(x) != length(y)) 
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stopO'Incorrect vector length") 
if (length(x) == 0) 
return(x) 

z <- match.fun(f )(x, y) 
return(aplRDV(z, g)) 

> 

# expand vector 

aplEXV <- function(x, y) { 
z <- rep(0, length(y)) 
m <- which (y == TRUE) 
if (length (m) != length(x)) 

stop ("Incorrect vector length") 
z[m] <- x 
return(z) 


# expand 

aplExpand <- function(x, y, axis = 1) { 
if (is.vector(x) ) 

return ( aplEXV (x, y)) 
d <- dim(x) 
m <- which (y == TRUE) 
n <- length (y) 
e <- d 

e[axis] <- n 

if (length (m) != d[axis]) 

stopO'Incorrect dimension length") 
z <- array (0, e) 
for (i in l:prod(d)) { 
k <- aplEncode (i, d) 
k[axis] <- m[k[axis]] 
z[aplDecode (k, e)] <- x[i] 

> 

return (z) 

> 

# compress/replicate vector 

aplCRV <- function(x, y) { 
n <- aplShape(x) 
m <- aplShape(y) 
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if (m == 1) 

y <- rep(y, n) 
if (length(y) != n) 
stopO'Length Error") 
z <- vectorO 
for (i in l:n) 

z <- c (z, rep(x [i] , y[i])) 
return(z) 


# compress/replicate 

aplReplicate <- function(x, y, k = aplRank (y)) { 
if (is.vector(x) ) 

return(aplCRV(x, y)) 
sx <- aplShape(x) 
sy <- aplShape(y) 
sk <- sx[k] 
if (max(sy) == 1) 
y <- rep(y , sk) 
if (length(y) != sk) 
stop ("Length Error") 
sz <- sx 
sz[k] <- sum(y) 
nz <- prod(sz) 
gg <- aplCRV(l :sk, y) 
z <- array(0, sz) 
for (i in l:nz) { 

jvec <- aplEncode(i, sz) 
jvec [k] <- gg [jvec [k] ] 
z [i] <- x [aplDecode(jvec, sx)] 

} 

return(z) 


# rotate vector 

aplRTV <- function(a, k) { 
n <- aplShape(a) 
if (k > 0) 

return(c(a[-(l :k)], a[l:k])) 
if (k < 0) 

return(c(a[(n + k + l):n], a[l:(n + k)])) 
return(a) 
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> 


# rotate 

aplRotate <- function(a, b, axis = aplRank (a)) { 
if (is.vector(a) ) 

return(aplRTV(a, b)) 
sa <- aplShape(a) 
sx <- aplShape(b) 
if (max(sx) == 1) 

b <- array(b, sa[-axis]) 
if (! identical (sa[-axis], aplShape (b))) 
stop ("Index Error") 
z <- array(0, sa) 
sz <- sa 
nz <- prod(sz) 
sk <- sz[axis] 
for (i in l:nz) { 

ivec <- aplEncode(i, sz) 
xx <- b [aplDecode (ivec[-axis], sx)] 
ak <- rep(0, sk) 
for (j in 1:sk) { 
jvec <- ivec 
jvec[axis] <- j 

ak[j] <- a [aplDecode (jvec, sz)] 

> 

bk <- aplRTV(ak, xx) 
for (j in 1:sk) { 
jvec <- ivec 
jvec[axis] <- j 

z [aplDecode (jvec, sz)] <- bk[j] 

> 

} 

return(z) 


# transpose — will be overwritten by the C version 

aplTranspose <- function(a, x = rev(l : aplRank(a) )) { 
sa <- aplShape (a) 
ra <- aplRank (a) 
if (length (x) != ra) 
stop ("Length Error") 
rz <- max(x) 
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sz <- rep(0, rz) 
for (i in l:rz) 

sz[i] <- min(sa[which(x == i)]) 
nz <- prod(sz) 
z <- array(0, sz) 
for (i in l:nz) 

z[i] <- a [aplDecode(aplEncode(i, sz) [x], sa)] 
return(z) 


# representation — will be overwritten by the C version 

aplEncode <- function(rrr, base) { 
b <- c(l, butLast(cumprod(base))) 
r <- rep(0, length(b)) 
s <- rrr - 1 

for (j in length(base) : 1) { 
r[j] <- s •/./•/. b[j] 
s <- s - r[j] * b[j] 

} 

return (1 + r) 


# base value — will be overwritten by the C version 

aplDecode <- function(ind, base) { 
b <- c(l, butLast(cumprod(base) )) 
return (1 + sum(b * (ind - 1))) 


# get 

aplGet <- function(a, cell) { 
dims <- dim(a) 
n <- length(dims) 
b <- 0 

if (any(cell > dims) || any(cell < 1)) 
stop("No such cell") 
return(a[aplDecode (cell, dims)]) 

> 

# set 

aplSet <- function(a, b, cell) { 
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dims <- dim(a) 
n <- length(dims) 

if (any(cell > dims) || any(cell < 1)) 
stopC'No such cell") 
a[aplDecode(cell, dims)] <- b 
return(a) 


# join 

aplJoin <- function(a, b, axis = 1) { 
if (is.vector(a) && is.vector(b) ) 
return(c(a, b)) 
sa <- aplShape(a) 
sb <- aplShape(b) 
ra <- aplRank(a) 
rb <- aplRank(b) 
if (ra != rb) 

stopO'Rank error in aplJoin") 
if (! identical(sa[-axis] , sb[-axis])) 
stop ("Shape error") 
sz <- sa 

sz[axis] <- sz[axis] + sbEaxis] 
nz <- prod(sz) 
u <- unit (axis, ra) 
z <- array(0, sz) 
for (i in l:nz) { 

ivec <- aplEncode(i, sz) 
if (ivecEaxis] <= saEaxis]) 

z[i] <- aEaplDecode(ivec, sa)] 
else 

z[i] <- bEaplDecode(ivec - saEaxis] * u, sb)] 

} 

return(z) 


# ravel 

aplRavel <- function(a) { 
as.vector(a) 

> 

# outer product 


33 


aplOuterProduct <- function(x, y, f = "*") { 
return(outer (x, y, f)) 

> 

# shape 

aplShape <- function(a) { 
if (is.vector(a) ) 
return(length(a) ) 
return(dim (a)) 


# rank 

aplRank <- function(a) { 
aplShape(aplShape (a)) 

> 

# reshape 

aplReshape <- function(a, d) { 
return(array(a, d)) 

> 

# reduce — will be overwritten by C version 

aplReduce <- function(a, k), f = "+") { 
if (is.vector(a) ) 

return(aplRDV(a, f)) 
ff <- if (is.function(f )) 
f 

else 

match.fun (f) 
sa <- aplShape (a) 
ra <- aplRank (a) 
na <- prod(sa) 
sz <- sa[(l:ra)[—k]] 
z <- array(0, sz) 
nz <- prod(sz) 
ind <- rep(0, nz) 
for (i in l:na) { 

ivec <- aplEncode(i, sa) 
jind <- aplDecode(ivec [—k], sz) 
if (ind[jind] == 0) { 
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zEjind] <- a[i] 
ind[jind] <- 1 

> 

else 

zEjind] <- ff (zEjind], aEi]) 

} 

return(z) 


# scan — will be overwritten by the C version 

aplScan <- function(a, k = aplRank(a) , f = "+") { 
if (is.vector(a) ) 

return(aplSCV(a, f)) 
ff <- if (is.function(f )) 
f 

else 

match.fun (f) 
sa <- aplShape(a) 
ra <- aplRank(a) 
sk <- saEk] 
u <- unit(k, ra) 
na <- prod(sa) 
z <- a 

for (i in l:na) { 

ivec <- aplEncode(i, sa) 
sk <- ivecEk] 
if (sk == 1) 
zEi] <- aEi] 
else 

zEi] <~ ff (z EaplDecode(ivec - u, sa)], aEi]) 

} 

return(z) 


# inner product — will be overwritten by the C version 

aplInnerProduct <- function(a, b, f = g = "+") { 

sa <- aplShape(a) 
sb <- aplShape(b) 
ra <- aplRank(a) 
rb <- aplRank(b) 
ia <- 1:(ra - 1) 
ib <- (ra - 1) + (l:(rb - 1)) 
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ff <- match.fun(f) 
gg <- match.fun(g) 
ns <- last(sa) 
nt <- first (sb) 
if (ns != nt) 

stop ("Incompatible array dimensions") 
sz <- c(butLast (sa), butFirst (sb)) 
nz <- prod(sz) 
z <- array(0, sz) 
for (i in l:nz) { 

ivec <- aplEncode(i, sz) 
for (j in l:ns) { 

aa <- a[aplDecode(c(ivec [ia], j), sa)] 
bb <- b [aplDecode(c(j , ivec[ib]), sb)] 
tt <- ff(aa, bb) 
if (j == 1) 
z[i] <- tt 
else 

z [i] <- gg(z [i] , tt) 

> 

} 

return(z) 


# member of 

aplMemberOf <- function(a, b) { 
sa <- aplShape(a) 
sb <- aplShape(b) 
na <- prod(sa) 
nb <- prod(sb) 
z <- array (0, sa) 
for (i in l:na) { 
aa <- a[i] 
for (j in l:nb) 
if (aa == b [j]) 
z[i] <- 1 

} 

return(z) 


# utilities below 
first <- function(x) { 
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return(x [1]) 


> 

butFirst <- function(x) { 
return(x [-1]) 

> 

last <- function(x) { 
return (x [length (x)]) 

> 

butLast <- function(x) { 
return(x[-length(x)] ) 

> 

unit <- functiond, n) { 
ifelse(i == (l:n), 1, 0) 

> 


5.2 R Glue for ,C() 


aplDecode <- 

function(cell, dims) { 
n <- length(dims) 

if (any (cell > dims) || any (cell < 1)) 
stopO'No such cell") 

.C( "aplDecodeC ", 
as.integer (cell) , 
as.integer (dims), 
as.integer (n), 
as.integer (1) )[ [4]] 


aplEncode <- 

functiondnd, dims) { 
n <- length(dims) 
cell <- integer(n) 

if ((ind < 1) || (ind > prod(dims))) 
stopC'No such cell") 

.C( "aplEncodeC ", 
as.integer (cell) , 
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as.integer (dims), 

as.integer (n), 

as.integer (ind))[[!]] 


aplInnerProduct <- 

function(a, b, f = g = "+") { 

sa <- aplShape(a) 
sb <- aplShape(b) 
ra <- aplRank(a) 
rb <- aplRank(b) 
ia <- 1 :(ra - 1) 
ib <- (ra - 1) + (l:(rb - 1)) 
ff <- match.fun(f) 
gg <- match.fun(g) 
ns <- last(sa) 
nt <- first(sb) 
if (ns != nt) 

stop ( "Incompatible array dimensions") 
sz <- c(butLast(sa) , butFirst (sb)) 
nz <- prod(sz) 
z <- array (0, sz) 
rz <- aplRank(z) 
res <- 

• C( 

"aplInnerProductC" , 
list (ff, gg), 
as.double (a), 
as.double (b), 
as.integer (sa), 
as.integer (ra), 
as.integer (sb), 
as.integer (rb), 
as.integer (sz), 
as.integer (rz), 
as.integer (nz), 
as.integer (ns), 
as.double (z) 

) 

return(array(res [ [12] ] , sz)) 

} 

aplReduce <- 
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function(a, k, f = "+") { 
if (is.vector(a) ) 

return(aplRDV (a, f)) 
ff <- if (is.function(f )) 
f 

else 

match.fun(f) 
sa <- aplShape(a) 
ra <- aplRank(a) 
na <- prod(sa) 
sz <- sa[ (1: ra) [—k] ] 
z <- array (0, sz) 
rz <- aplRank(z) 
nz <- prod(sz) 
z <- 
• C( 

"aplReduceC" , 

list (ff) , 

as.double (a), 

as.integer (k), 

as.integer(length(k) ), 

as.integer (na), 

as.integer (sa), 

as.integer (ra), 

as.integer (nz), 

as.integer (sz), 

as.integer (rz), 

as.double (z) 

) 

return(array(z [ [11] ], sz)) 

} 

aplScan <- 

function(a, k, f = "+") { 
if (is.vector (a)) 

return(aplSCV(a, f)) 
ff <- if (is.function(f )) 
f 

else 

match.fun (f) 
sa <- aplShape(a) 
ra <- aplRank(a) 
sk <- sa[k] 
u <- unit(k, ra) 
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na <- prod(sa) 
z <- a 
res <- .C( 

"aplScanC" , 
list (ff), 
as.double (a), 
as.integer (k), 
as.integer (na), 
as.integer (sa), 
as.integer (ra), 
as.double (z) 

) 

return(array(res [[7]], sa)) 

} 

aplSelect <- 

function(a, x, drop = FALSE) { 
sa <- aplShape(a) 
ra <- aplRank(a) 
sz <- sapply(x, length) 
z <- array (0, sz) 
rz <- aplRank(z) 
nz <- prod(sz) 
z <- 
array ( 

• C( 

"aplSelectC" , 

as.double (a), 

as.integer (sa), 

as.integer (ra), 

lapply(x, as.integer), 

as.double (z), 

as.integer (sz), 

as.integer (rz), 

as.integer (nz) 

)[[5]] , 

sz 

) 

if (drop) 

return(drop(z) ) 
else 

return (z) 
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aplTranspose <- 

function(a, x = rev(l : aplRank(a) )) { 
sa <- aplShape(a) 
ra <- aplRank(a) 
na <- prod(sa) 
if (length(x) != ra) 
stop ("Length Error") 
rz <- max(x) 
sz <- rep(0, rz) 
for (i in l:rz) 

sz[i] <- min(sa[which(x == i)]) 
nz <- prod(sz) 
z <- array (0, sz) 
array ( 

• C( 

"aplTransposeC" , 
as.double (a), 
as.integer (x), 
as.integer (sa), 
as.integer (ra), 
as.integer (na), 
as.integer (sz), 
as.integer (rz), 
as.integer (nz), 
as.double (z) 

) [[9]] , 

sz 

) 

} 

5.3 R Glue for .CallO 


aplDecode <- function(cell, dims) { 
if (length(cell) != length(dims) ) { 
stopC'Dimension error") 

> 

if (any(cell > dims) || any (cell < 1)) { 
stopO'No such cell") 

> 

.Call ( "APLDECODE" , as.integer (cell), as.integer (dims)) 
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aplEncode <- function(ind, dims) { 
if (length(ind) > 1) { 
stop ("Dimension error") 

> 

if ((ind < 1) || (ind > prod(dims))) { 
stop ("No such cell") 

> 

.Call ( "APLENCODE" , as.integer (ind), as.integer (dims)) 


aplInnerProduct <- function(a, b, f = g = "+") { 

sa <- aplShape(a) 
ra <- aplRank(a) 
ia <- 1:(ra - 1) 
sb <- 

aplShape(b) 
rb <- aplRank(b) 
ib <- (ra - 1) + (l:(rb - 1)) 
if (! is.function(f )) { 
f <- match.fun(f) 

} 

if (! is.function(g) ) { 
g <- match.fun(g) 

} 

env <- new.envO 
environment (f) <- env 
environment (g) <- env 
ns <- last(sa) 
nt <- first (sb) 
if (ns != nt) { 

stop ("Incompatible array dimensions") 

} 

sz <- c(butLast (sa), butFirst (sb)) 
z <- . Call ( 

"APLINNERPRODUCT", 

f, 
g> 

as.double (a), 
as.double (b), 
as.integer (sa), 
as.integer (sb), 
as.integer (sz), 
as.integer (ns), 
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env 


) 

return(array (z, sz)) 


aplReduce <- function(a, k, f = "+") { 
if (is.vector(a) ) { 
return(aplRDV(a, f)) 

} 

nk <- length(k) 
sa <- aplShape(a) 
sz <- sa[(1: length(sa) )[—k]] 
if (! is.function(f )) { 
f <- match.fun(f) 

} 

env <- new.env () 
environment (f) <- env 
z <- .Call ( "APLREDUCE" , 

f, 

as.double (a), 
as.integer (k), 
as.integer (sa), 
as.integer (sz), 
env) 

return(array (z, sz)) 

> 

aplScan <- function(a, k = aplRank (a), f = "+") { 
if (is.vector(a) ) { 
return(aplSCV(a, f)) 

} 

if (! is.function(f )) { 
f <- match.fun(f) 

} 

env <- new.env () 
environment (f) <- env 
sa <- aplShape(a) 
ra <- aplRank (a) 
z <- .Call ( "APLSCAN" , 

f, 

as.double (a), 
as.integer (k), 
as.integer (sa), 
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as.integer (ra), 
env) 

return(array (z, sa)) 

> 

aplSelect <- function (a, x, drop = TRUE) { 
dima = aplShape(a) 
if (length(dima) != length(x)) { 
stopC'Dimension error") 

} 

z <- .Call ( "APLSELECT" , 
as.double (a), 
as.integer (dima), 
lapply(x, as.integer)) 

z <- array(z, sapply(x, length)) 
if (drop) { 

return(drop(z) ) 

} 

return(z) 


aplTranspose <- function(a, x = rev(l : aplRank(a) )) { 
sa <- aplShape(a) 
if (length(x) != length(sa)) { 
stopC'Length Error") 

} 

rz <- max(x) 
sz <- rep(0, rz) 
for (i in l:rz) { 

sz[i] <- min(sa[which(x == i)]) 

} 

z <- .Call( 

"APLTRANSPOSE", 
as.double (a), 
as.integer (x), 
as.integer (sa), 
as.integer (sz), 
as.integer (rz) 

) 

return(array (z, sz)) 
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5.4 C Code using ,C() 


#include <R.h> 

#include <Rinternals,h> 

#include <Rinterface,h> 

static char* f2_char; 
static char* g2_char; 

void aplDecodeC(int* cell, int* dims, int* n, int* ind) { 
int i, aux = 1; *ind = 1; 
for (i = 0; i < *n; i++) { 

*ind += aux * (cell[i] - 1); 
aux *= dims[i]; 

} 

} 

void aplEncodeC(int* cell, int* dims, int* n, int* ind) { 
int i, aux = *ind, nn = *n, pdim = 1; 
for (i = 0; i < nn - 1; i++) 
pdim *= dims[i]; 
for (i = nn - 1; i > 0; i—) { 
cell [i] = (aux - l)/pdim; 
aux -= pdim*cell[i]; 
pdim /= dims[i-l]; 
cell [i] += 1; 

} 

cell[0] = aux; 

} 

void aplSelectC (double *a, int *sa, int *ra, SEXP *list, double *z, int *sz, int *rz, in 
int i, j, k; 

int ivec[*rz], jvec[*ra]; 
for (i = 0; i < *nz; i++) { 
k = i + 1; 

(void) aplEncodeC (ivec,sz,rz,&k); 
for (j = 0; j < *ra; j++) { 

SEXP vec = list[j]; 

jvec[j] = INTEGER(vec) [ivec [j] -1] ; 

} 

k = 1; 

(void) aplDecodeC (jvec,sa,ra,&k); 
z [i] = a[k - 1] ; 
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} 


> 

void aplTransposeC (double *a, int *x, int *sa, int *ra, int *na, int *sz, int *rz, int * 

int i, j, r, ivec [*rz], jvec[*ra]; 
for (i = 0; i < *nz; i++){ 
r = i + 1; 

(void) aplEncodeC(ivec,sz,rz,&r); 
for (j = 0; j < *ra; j++) 
jvec [j] = ivec [x [j] -1] ; 

(void) aplDecodeC(jvec,sa,ra,&r); 
z [i] = a[r - 1] ; 

} 


void aplReduceC(char** funclist, double *a, int *k, int *nk, int *na, int *sa, int *ra, 
double f2_glue (double, double); 
double f2_comp (double (*)(), double,double) ; 

int i, j, u, v, r, m, kk= (*ra) - (*nk), ivec[*ra], kvec[kk], ind[*nz]; 
for (i = 0; i < *nz; i++) ind[i] = 0; 
f2_char = funclist[0]; 
for (i = 0; i < *na; i++){ 
r = i + 1; 

(void) aplEncodeC(ivec,sa,ra,&r); 
u = 0; 

for (j = 0; j < *ra; j++) { 
r = 0; 

for (v = 0; v < *nk; v++) { 
if (j == (k[v] - 1)) r = 1; 

} 

if (r == 0) { 

kvec[u] = ivec[j]; 
u += 1; 

} 

> 

(void) aplDecodeC(kvec,sz,rz,&m); 
if (ind[m - 1] == 0) { 
z [m - 1] = a[i] ; 
ind[m - 1] =1; 

} 

else 

z[m - 1] = f2_comp( (double (*)())f2_glue,z[m - 1] ,a[i]); 

} 
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> 


void aplScanC(char** funclist, double *a, int *k, int *na, int *sa, int *ra, double *z) 

double f2_glue (double, double); 

double f2_comp (double (*)(), double,double) ; 

int i, r, sk, 1, ivec[*ra]; 

1 = *k - 1; 

f2_char = funclist[0]; 
for (i = 0; i < *na; i++){ 
r = i + 1; 

(void) aplEncodeC(ivec,sa,ra,&r); 
sk = ivec [1]; 
if (sk == 1) z[i] = a[i] ; 
else { 

ivec [1] -= 1; 

(void) aplDecodeC(ivec,sa,ra,&r); 

z[i] = f2_comp( (double (*)())f2_glue,z[r - 1] ,a[i]); 

> 

} 

> 

void aplInnerProductC(char** funclist, double *a, double *b, int *sa, int *ra, int *sb, 
double f2_glue (double, double); 
double g2_glue (double, double); 
double f2_comp (double (*)(), double,double) ; 
double g2_comp (double (*)(), double,double) ; 

double t; int i, j, r, k, 1, u, ivec[*rz], jvec[*ra], kvec[*rb]; 
f2_char = funclist[0]; g2_char= funclist[1]; k = 1 = 0; 
for (i = 0; i < *nz; i++) { 
r = i + 1; 

(void) aplEncodeC(ivec,sz,rz,&r); 
for (j = 0; j < *ns; j++) { 
for (u = 0; u < *ra - 1; u++) 
jvec[u] = ivec [u] ; 
jvec[*ra - 1] = j + 1; 

(void) aplDecodeC(jvec,sa,ra,&k); 
for (u = 1; u < *rb; u++) 
kvec[u] = ivec[*ra + u - 2] ; 
kvec[0] = j + 1; 

(void) aplDecodeC(kvec,sb,rb,&l); 
t = f2_comp( (double (*) 0)f2_glue,a[k-l] ,b[1—1]) ; 
if (j == 0) z [i] = t; 

else z[i] = g2_comp( (double (*) Q)g2_glue,t,z[i]) ; 
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> 


} 

> 

double f2_glue (double x, double y){ 
char *modes [2] , *values[l]; 
void *args[2]; 

double xx[1], yy[l], *result; 

long lengths [2], nargs = (long) 2, nvals = (long)l; 
lengths[0] = lengths[1] = (long)l; 
nargs = (long) 2; 

args[0] = (void *)xx; xx[0] = x; 
args[1] = (void *)yy; yy[0] = y; 
modes[0] = modes[1] = "double"; 

call_R(f2_char,nargs,args,modes,lengths, 0 ,nvals,values); 
result = (double*) values[0]; 
return(result [0]); 

> 

double g2_glue (double x, double y){ 
char *modes [2] , *values[l]; 
void *args[2]; 

double xx [1], yy[l], *result; 

long lengths[2], nargs = (long) 2, nvals = (long)l; 
lengths[0] = lengths[1] = (long)l; 
nargs = (long) 2; 

args[0] = (void *)xx; xx[0] = x; 
args[1] = (void *)yy; yy[0] = y; 
modes[0] = modes[1] = "double"; 

call_R(g2_char,nargs,args,modes,lengths, 0 ,nvals,values); 
result = (double*) values[0]; 
return(result [0]); 

> 

double f2_comp (double (*f)(), double x, double y) { 
return f(x,y); 

> 

double g2_comp (double (*g)(), double x, double y) { 
return g(x,y); 

> 


48 


5.5 C Code using .Call() 


#include <R.h> 

#include <Rinternals,h> 

SEXP APLDECODE( SEXP, SEXP ); 

SEXP APLENCODE( SEXP, SEXP ); 

SEXP APLSELECT( SEXP, SEXP, SEXP ); 

SEXP APLTRANSPOSE( SEXP, SEXP, SEXP, SEXP, SEXP ); 

SEXP APLSCAN( SEXP, SEXP, SEXP, SEXP, SEXP ); 

SEXP APLREDUCE( SEXP, SEXP, SEXP, SEXP, SEXP, SEXP ); 

SEXP APLINNERPRODUCT( SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP ); 

SEXP 

APLDECODE( SEXP cell, SEXP dims ) 

{ 

int aux = 1, n = length(dims); 

SEXP ind; 

PROTECT( ind = allocVector( INTSXP, 1 ) ); 

INTEGER( ind )[0] = 1; 

for( int i = 0; i < n; i++ ) { 

INTEGERC ind )[0] += aux * ( INTEGER( cell )[i] - 1 ); 
aux *= INTEGER(dims)[i]; 

> 

UNPROTECT( 1 ); 
return (ind); 

} 

SEXP 

APLENCODE( SEXP ind, SEXP dims ) 

{ 

int n = length(dims), aux = INTEGER(ind)[0], pdim = 1; 

SEXP cell; 

PROTECT( cell = allocVector( INTSXP, n ) ); 
for ( int i=0; i < n - 1; i++ ) 
pdim *= INTEGER( dims )[i]; 
for ( int i=n-l;i>0; i— ){ 

INTEGER( cell )[i] = ( aux - 1 ) / pdim; 
aux -= pdim * INTEGER( cell )[i]; 
pdim /= INTEGER( dims )[i - 1]; 

INTEGER( cell )[i] += 1; 

> 

INTEGER( cell )[0] = aux; 

UNPROTECT( 1 ); 
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> 


return cell; 


SEXP 

APLSELECTC SEXP a, SEXP dima, SEXP list ) { 

int r = length (dima), lz = 1 , dimzi, nProtect = 0; 


SEXP 

PROTECT(dimz 
PROTECT(cell 
PROTECT(czll 
PROTECT(itel 
PROTECT(nind 
for( int i = 


dimz, itel, cell, czll, 

= allocVector (INTSXP, r)) 
= allocVector (INTSXP, r)) 
= allocVector (INTSXP, r)) 
= allocVector (INTSXP, 1)) 
= allocVector (INTSXP, 1)) 
0; i < r; i++ ){ 


nind, z; 

; nProtect++ 
; nProtect++ 
; nProtect++ 
; nProtect++ 
; nProtect++ 


dimzi = length( VECT0R_ELT( list, i ) ); 
INTEGER( dimz )[i] = dimzi; 
lz *= dimzi; 


> 

PROTECT( z = allocVector( REALSXP, lz ) ); 
nProtect++; 

for( int i = 0; i < lz; i++ ){ 

INTEGER(itel)[0] = i + 1; 
cell = APLENCODE (itel, dimz); 
for (int j = 0; j < r; j++) { 

INTEGER (czll) [j] = INTEGER (VECT0R_ELT (list, j))[INTEGER(cell)[j] - 1]; 

} 

nind = APLDEC0DE( czll, dima ); 

REAL(z) [i] = REAL(a)[INTEGER(nind)[0] - 1]; 


> 

UNPROTECT(nProtect); 
return(z) ; 


} 


SEXP 

APLTRANSP0SE(SEXP a, SEXP x, SEXP sa, SEXP sz, SEXP rz) 

{ 


int i, j, na=l, nz=l, ra = length (sa), lsz = length( sz 

SEXP ivec, jvec, z, itel, nind; 

for( i=0;i<ra ;i++){ na *= INTEGER( sa ) [i]; } 

for( i=0;i<lsz;i++){ nz *= INTEGER( sz ) [i]; } 

PROTECT(itel = allocVector( INTSXP, 1 ) ); 

PROTECT(nind = allocVector( INTSXP, 1 ) ); 

PROTECT(ivec = allocVector( INTSXP, INTEGER(rz)[0] ) ); 
PROTECT(jvec = allocVector( INTSXP, ra ) ); 

PROTECT(z = allocVector( REALSXP, nz ) ); 


), nProtected=0; 


++nProtected; 

++nProtected; 

++nProtected; 

++nProtected; 

++nProtected; 
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for( i = 0; i < nz; i++ ){ 

INTEGER( itel )[0] = i + 1; 
ivec = APLENC0DE( itel, sz ); 
for ( j = 0; j < ra; j++ ){ 

INTEGERC jvec )[j] = INTEGER( ivec ) [INTEGER( x )[j] - 1] ; 

} 

nind = APLDEC0DE( jvec, sa ); 

REAL( z )[i] = REAL( a )[INTEGER(nind) [0] - 1] ; 

> 

UNPROTECT( nProtected ); 
return z; 


SEXP 

APLREDUCE(SEXP f, SEXP a, SEXP k, SEXP sa, SEXP sz, SEXP env) 

{ 

int i, j, u, v, r, kk, na=l, nz=l, nProtected = 0; 

int nk = length( k ), ra = length( sa ), rz = length( sz ); 

SEXP ivec, kvec, ind, z,itel, nind, R_fcall= R_NilValue; 

SEXP Z = R_NilValue, A = R_NilValue; 

for( i=0;i<ra;i++){ na *= INTEGERC sa ) [i]; } 

for( i=0;i<rz;i++){ nz *= INTEGERC sz ) [i]; } 

kk = ra-nk; 


PROTECT( 

R_fcall= lang3(f, R_NilValue, 

R_NilValue) 

); 

++nProtected 

PROTECT( 

Z = allocVectorC 

REALSXP, 

1 

) 

); 

++nProtected 

PROTECT( 

A = allocVectorC 

REALSXP, 

1 

) 

); 

++nProtected 

PROTECT( 

itel = allocVectorC 

INTSXP, 

1 

) 

); 

++nProtected 

PROTECT( 

ivec = allocVectorC 

INTSXP, 

ra 

) 

); 

++nProtected 

PROTECT( 

kvec = allocVectorC 

INTSXP, 

kk 

) 

); 

++nProtected 

PROTECT( 

ind = allocVectorC 

INTSXP, 

nz 

) 

); 

++nProtected 

PROTECT( 

z = allocVectorC 

REALSXP, 

nz 

) 

); 

++nProtected 


for( i = 0; i < nz; i++ ){ 

INTEGERC ind )[i] = 0; 

> 

for( i = 0; i < na; i++ ){ 

INTEGERC itel )[0] = i + 1; 
ivec = APLENCODEC itel, sa ); 
u = 0; 

for( j = 0; j < ra; j++ ){ 
r = 0; 

for( v = 0; v < nk; v++ ){ 

if ( j == ( INTEGERC k )[v] - 1 ) ) r = 1; 

} 

if ( r == 0 ){ 
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INTEGERS kvec )[u] = INTEGER( ivec )[j]; 
u += 1 ; 

} 

} 

nind = APLDECODE( kvec, sz ); 

if ( INTEGER( ind )[INTEGER( nind )[0] - 1] == 0 ) { 

REAL( z )[INTEGER( nind )[0] - 1] = REAL( a )[i]; 

INTEGER( ind )[INTEGER( nind )[0] - 1] =1; 

} else { 

REAL( Z )[0] = REAL( z )[INTEGER( nind )[0] - 1]; 

REAL( A ) [0] = REAL( a ) [i] ; 

SETCADRC R_fcall, Z ); 

SETCADDRC R_fcall, A ); 

REAL( z )[INTEGER( nind )[0] - 1] = REAL( eval( R_fcall, env ) )[0]; 

} 

> 

UNPROTECT( nProtected ); 
return z; 

> 


SEXP 

APLSCAN( SEXP f, SEXP a, SEXP k, SEXP sa, SEXP ra, SEXP env ) 

{ 


int i, sk, 1, na=l, ka = INTEGER( ra )[0],nProtected=0; 
for( i=0;i<ka ;i++ ){ na *= INTEGER( sa ) [i]; } 

SEXP ivec, z, itel, nind, Z, A,R_fcall=R_NilValue; 
PROTECT( R_fcall = lang3( f, R_NilValue, R_NilValue ) ) 


PROTECT( itel = allocVector( INTSXP, 1 ) ) 
PROTECT( nind = allocVector( INTSXP, 1 ) ) 
PROTECT( Z = allocVector( REALSXP, 1 ) ) 
PROTECT( A = allocVector( REALSXP, 1 ) ) 
PROTECT( ivec = allocVector( INTSXP, ka ) ) 
PROTECT( z = allocVector( REALSXP, na ) ) 


1 = INTEGERC k )[0] - 1; 
for( i = 0; i < na; i++ ){ 
INTEGERC itel )[0] = i + 1; 
ivec = APLENCODEC itel, sa ); 
sk = INTEGERC ivec )[1]; 
if ( sk == 1 ){ 


++nProtected; 

++nProtected; 

++nProtected; 

++nProtected; 

++nProtected; 

++nProtected; 

++nProtected; 


REALC z ) [i] = REALC a ) [i] ; 

}else{ 

INTEGERC ivec )[1] -= 1; 
nind = APLDEC0DEC ivec, sa ); 

REALC Z ) [0]=REAL( z )[INTEGERC nind ) [0] -1] ; 
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REAL( A )[0]=REAL( a )[i]; 

SETCADRC R_fcall, Z ); 

SETCADDRC R_fcall, A ); 

REAL( z )[i] = REAL( eval( R_fcall, env ) ) [0]; 

} 

> 

UNPROTECT( nProtected ); 
return z; 


SEXP 

APLINNERPRODUCT(SEXP f, SEXP g, SEXP a, SEXP b, SEXP sa, SEXP sb, SEXP sz, SEXP ns, SEXF 
int i, j, u, nz = 1, nProtected = 0; 

int ra = length( sa ),rb = length( sb ), rz = length( sz ); 

SEXP ivec, jvec, kvec, z, A, B, Z, itel, k, 1, t; 

SEXP R_fcall = R_NilValue, R_gcall = R_NilValue; 


forC i = 

0; i < 

rz; i++ ){ nz 

*= INTEGERC 

sz ) [i] ; 

> 



PROTECT(R_fcall 

= 

lang3( f, R_ 

lilValue, 

R. 

NilValue 

) 

) 

++nProtected 

PROTECT(R_gcall 

= 

lang3( g, R_ 

lilValue, 

R. 

NilValue 

) 

) 

++nProtected 

PROTECT( 

itel 

= 

allocVector( 

INTSXP, 


1 

) 

) 

++nProtected 

PROTECT( 

k 

= 

allocVector( 

INTSXP, 


1 

) 

) 

++nProtected 

PROTECTC 

1 

= 

allocVectorC 

INTSXP, 


1 

) 

) 

++nProtected 

PROTECTC 

ivec 

= 

allocVectorC 

INTSXP, 


rz 

) 

) 

++nProtected 

PROTECTC 

jvec 

= 

allocVectorC 

INTSXP, 


ra 

) 

) 

++nProtected 

PROTECTC 

kvec 

= 

allocVectorC 

INTSXP, 


rb 

) 

) 

++nProtected 

PROTECTC 

z 

= 

allocVector( 

REALSXP, 


nz 

) 

) 

++nProtected 

PROTECT( 

t 

= 

allocVector( 

REALSXP, 


1 

) 

) 

++nProtected 

PROTECT( 

Z 

= 

allocVector( 

REALSXP, 


1 

) 

) 

++nProtected 

PROTECT( 

B 

= 

allocVector( 

REALSXP, 


1 

) 

) 

++nProtected 

PROTECT( 

A 

= 

allocVector( 

REALSXP, 


1 

) 

) 

++nProtected 


for( i = 0; i < nz; i++ ){ 

INTEGERC itel )[0] = i + 1; 
ivec = APLENCODE( itel, sz ); 
for( j = 0; j < INTEGERC ns ) [0]; j++ ){ 
for( u=0;u<ra-l; u++ ){ 

INTEGERC jvec )[u] = INTEGERC ivec )[u]; 

} 

INTEGERC jvec )[ra - 1] = j + 1; 
k = APLDECODEC jvec, sa ); 
for( u = 1; u < rb; u++ ){ 

INTEGERC kvec )[u] = INTEGERC ivec )[ra + u - 2] ; 

} 

INTEGERC kvec )[0] = j + 1; 
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1 = APLDECODE( kvec,sb ); 

REAL( A )[0] = REAL( a )[INTEGER( k ) [0] -1] ; 

REAL( B )[0] = REAL( b )[INTEGER( 1 ) [0] -1] ; 
SETCADRC R_fcall, A ); 

SETCADDRC R_fcall, B ); 

REAL( t )[0] = REAL( eval( R_fcall, env ) ) [0] ; 
if ( j == 0 ){ 

REAL( z ) [i] = REAL( t ) [0] ; 

} else { 

REAL( Z ) [0] = REAL( z ) [i] ; 

SETCADRC R_gcall, t ); 

SETCADDRC R_gcall, Z ); 

REAL( z )[i] = REAL( evalC R_fcall, env ) ) [0]; 

} 

} 

> 

UNPROTECT( nProtected ); 
return z; 

> 


5.6 Repp code 
5.6.1 hpp 


#ifndef _apl_RCPP_H 
#define _apl_RCPP_H 

/* The Repp header takes care of everything you should need. */ 

#include <Rcpp.h> 

/* 

* note : RcppExport is an alias to 'extern "C"' defined by Repp. 

* 

* It gives C calling convention to the rcpp_hello_world function so that 

* it can be called from .Call in R. Otherwise, the C++ compiler mangles the 

* name of the function and .Call can't find it. 

* 

* It is only useful to use RcppExport when the function is intended to be called 

* by .Call. See the thread http: //thread, gmane. org/gmane. comp. lang. r. rcpp/6f.9/focus=6 

* on Rcpp-devel for a misuse of RcppExport 
*/ 

//RcppExport SEXP rcpp_hello_world() ; 
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RcppExport 

RcppExport 

RcppExport 

RcppExport 

RcppExport 

RcppExport 

RcppExport 


SEXP APLDECODErcpp( SEXP, SEXP ); 

SEXP APLENCODErcpp( SEXP, SEXP ); 

SEXP APLSELECTrcpp( SEXP, SEXP, SEXP ); 

SEXP APLTRANSPOSErcpp( SEXP, SEXP, SEXP, SEXP, SEXP ); 

SEXP APLSCANrcppC SEXP, SEXP, SEXP, SEXP, SEXP ); 

SEXP APLREDUCErcpp( SEXP, SEXP, SEXP, SEXP, SEXP, SEXP ); 

SEXP APLINNERPRODUCTrcpp( SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEX 


#endif 


5.6.2 cpp 


// 

#include "aplRcpp.hpp" 
using namespace Repp ; 

SEXP 

APLDECODErcpp( SEXP celljr, SEXP dims_r ) 

{ 

IntegerVector cell( cell_r ); 

IntegerVector dims( dims_r ); 
int n = dims.size(), aux = 1; 

IntegerVector ind( 1 ); ind[0] = 1; 
for( int i = 0; i < n; i++ ) { 

ind[0] += aux * ( cell [i] - 1 ); 
aux *= dims[i]; 

} 

return wrap( ind ); 

} 

SEXP 

APLENCODErcpp( SEXP ind_r, SEXP dims_r ) 

{ 

IntegerVector dims( dims_r ); 

IntegerVector ind( ind_r ); 
int n = dims.size(), aux = ind[0], pdim = 1; 
IntegerVector cell( n ); 
for( int i=0;i<n-l; i++ ){ 
pdim *= dims[i]; 

> 

for( int i = n - 1; i>0; i— ){ 

cell [i] = ( aux - 1 ) / pdim; 
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aux -= pdim * cell [i]; 
pdim /= dims[i - 1]; 
cell [i] += 1; 

> 

cell [0] = aux; 
return wrap( cell ); 


SEXP 

APLSELECTrcpp( SEXP a_r, SEXP dima_r, SEXP list_r ) 

{ 


IntegerVector dima( dima_r ); 

NumericVector a( a_r ); 

int r = dima.sizeO, lz = 1, dimzi; 

List list( list_r ); 

IntegerVector dimz( r ); 

IntegerVector cell( r ); 

IntegerVector czll( r ); 

IntegerVector itel( 1 ); 

IntegerVector nind( 1 ); 
for( int i = 0; i < r; i++ ){ 

dimzi = as<IntegerVector>( list[i] ).size(); 
dimz[i] = dimzi; 


lz *= dimzi; 

> 

NumericVector z(lz); 
for( int i = 0; i < lz; i++ ){ 
itel [0] = i + 1 ; 

cell = APLENCODErcppC wrap( itel ), wrap( dimz ) ); 
for ( int j =0; j <r; j++) { 

IntegerVector listj= as<IntegerVector>( list[j] 
czll[j] = list j [cell [j] - 1]; 

} 

nind = APLDECODErcppC wrap( czll ), wrap( dima ) ); 
z [i] = a [nind [0] - 1]; 

> 

return wrap( z ); 


); 


SEXP 

APLTRANSPOSErcpp( SEXP a_r, SEXP x_r, SEXP sa_r, SEXP sz_r, SEXP rz_r ) 

IntegerVector sa( sa_r ), sz( sz_r ); 

int na = 1, nz = 1, ra = sa.sizeQ, lsz = sz.sizeO; 
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IntegerVector rz( rz_r ), x( x_r ), a( a_r ); 
for( int i=0;i< ra - 1; i++ ){ na *= sa[i]; } 

for( int i = 0; i < lsz - 1; i++ ){ nz *= sz[i]; } 

IntegerVector itel( 1 ), nind( 1 ), ivec( rz [0] ), jvec( ra ); 

NumericVector z( nz ); 
for( int i = 0; i < nz; i++ ){ 
itel[0] = i + 1; 

ivec = APLENCODErcppC wrap( itel ), wrap( sz ) ); 
for( int j=0;j<ra;j++){ 
jvec [j] = ivec [x [j] - 1] ; 

} 

nind = APLDECODErcppC wrap( jvec ), wrap( sa ) ); 
z[i] = a [nind [0] - 1]; 

> 

return wrap( z ); 

} 

SEXP 

APLREDUCErcpp( SEXP f_r, SEXP a_r, SEXP k_r, SEXP sa_r, SEXP sz_r, SEXP env_r ) 

{ 

int u, r, kk, na = 1, nz = 1, nProtected=0; 

IntegerVector k( k_r ), sa( sa_r ), sz( sz_r ); 

NumericVector a( a_r ); 

int nk = k.sizeO, ra = sa.sizeO, rz = sz.sizeO; 

SEXP R_fcall= R_NilValue; 

for( int i = 0; i < ra; i++ ){ na *= sa[i]; } 
for( int i = 0; i < rz; i++ ){ nz *= sz [i]; } 
kk = ra - nk; 

PR0TECT( R_fcall= Rf_lang3(f_r, R_NilValue, R_NilValue) ); ++nProtected; 
IntegerVector itel( 1 ), nind( 1 ), ind( nz ), ivec( ra ), kvec( kk ); 
NumericVector z( nz ); 
for( int i = 0; i < na; i++ ){ 
itel[0] = i + 1; 

ivec = APLENCODErcppC wrap( itel ), wrap( sa ) ); 
u = 0; 

for( int j = 0;j<ra;j++)-[ 
r = 0; 

for( int v = 0; v < nk; v++ ){ 

if( j == ( k[v] - 1 ) ) r = 1; 

} 

if( r == 0 ){ 

kvec[u] = ivec[j]; 
u += 1 ; 

} 
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} 

nind = APLDECODErcppC wrap( kvec ), wrap( sz ) ); 
if ( ind[nind[0] - 1] == 0 ) { 
z[nind[0] - 1] = a [i] ; 
ind[nind[0] - 1] = 1; 

} else { 

SETCADRC R_fcall, wrap( z[nind[0] - 1] ) ); 

SETCADDRC R_fcall, wrap( a[i] ) ); 

z[nind[0] - 1] = REAL( Rf_eval( R_fcall, env_r ) )[0] 

} 

> 

UNPROTECT( nProtected ); 
return wrap( z ); 

} 

SEXP 

APLSCANrcpp( SEXP f_r, SEXP a_r, SEXP k_r, SEXP sa_r, SEXP env_r 

{ 

IntegerVector k( k_r ) , sa( sa_r ); 

NumericVector a( a_r ); 

int sk, 1, na=l, ra = sa.sizeO, nProtected=0; 
for( int i = 0; i < ra; i++ ){ na *= sa[i]; } 

SEXP R_fcall=R_NilValue; 

PROTECT( R_fcall = Rf_lang3( f_r, R_NilValue, R_NilValue ) ); 
IntegerVector itel( 1 ), nind( 1 ), ivec( ra ); 

NumericVector z( na ); 

1 = k [0] - 1; 

for( int i = 0; i < na; i++ ){ 
itel[0] = i + 1; 

ivec = APLENCODErcppC itel, sa ); 
sk = ivec [1]; 
if( sk == 1 ){ 

z [i] = a[i] ; 

} else { 

ivec [1] -= 1; 

nind = APLDECODErcppC wrap( ivec ), wrap( sa ) ); 
SETCADRC R_fcall, wrap( z[nind[0]-l] ) ); 

SETCADDRC R_fcall, wrap( a[i] ) ); 

z[i] = REAL( Rf_eval( R_fcall, env_r ) )[0]; 

} 

> 

UNPROTECT( nProtected ); 
return wrap( z ); 

> 


++nProtected; 
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SEXP 

APLINNERPRODUCTrcpp(SEXP f_r, SEXP g_r, SEXP a_r, SEXP b_r, SEXP sa_r, SEXP sb_r, SEXP 

{ 


int nz = 1, nProtected = 0; 

IntegerVector sa( sa_r ) , sb( sb_r ) , sz( sz_r ), ns( ns_r ); 

NumericVector a( a_r ), b( b_r ); 

int ra = sa.size(),rb = sb.sizeO, rz = sz.sizeO; 

SEXP R_fcall = R_NilValue, R_gcall = R_NilValue; 
for( int i = 0; i < rz; i++ ){ nz *= sz [i]; } 

PR0TECT( R_fcall = Rf_lang3( f_r, R_NilValue, R_NilValue ) ); ++nProtected; 
PR0TECT( R_gcall = Rf_lang3( g_r, R_NilValue, R_NilValue ) ); ++nProtected; 
IntegerVector itel( 1 ), k( 1 ), 1( 1 ), ivec( rz ), jvec( ra ), kvec( rb ); 

NumericVector z( nz ), t( 1 ), Z( 1 ), B( 1 ), A( 1 ); 

for( int i = 0; i < nz; i++ ){ 
itel[0] = i + 1; 

ivec = APLENCODErcppC itel, sz ); 
for( int j = 0; j < ns [0]; j++ ){ 

for( int u=0; u < ra - 1; u++ ){ 
jvec[u] = ivec [u] ; 

} 

jvec[ra - 1] = j + 1; 
k = APLDECODErcppC jvec, sa ); 
for( int u = 1; u < rb; u++ ){ 
kvec [u] = ivec[ra + u - 2] ; 

} 

kvec [0] = j + 1; 

1 = APLDECODErcppC kvec,sb ); 

SETCADRC R_fcall, wrap( a[k[0]-l]) ); 

SETCADDRC R_fcall, wrap( b[1 [0] -1] )); 
t[0] = REAL( Rf_eval( R_fcall, env_r ) )[0]; 
if ( j == 0 ){ 
z [i] = t [0] ; 

} else { 

SETCADRC R_gcall, wrap( t[0] ) ); 

SETCADDRC R_gcall, wrap( z[i] ) ); 
z[i] = REAL( Rf_evalC R_gcall, env_r ) )[0]; 

} 

} 

> 

UNPROTECT( nProtected ); 
return wrap( z ); 

> 
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6 NEWS 


001 03/04/16 

• First release 
002 03/04/16 

• Make pdf, using Code2000 Unicode font 

• fixed bug in aplDropO 

003 03/05/16 

• fixed bug in aplExpandO 

• examples and text drop, expand, decode, encode 

• aplExpandO can now have logical or binary second argument 

004 03/05/16 

• bug corrected in APLINNERPRODUCT (.Call() interface). 

• writeup and examples aplExpandO 

• writeup and examples aplInnerProductO 

• bug corrected in aplJoinO 

• writeup and examples aplJoinO 

• writeup and example aplMemberOfO 

• changed aplMemberOfO so it applies to vectors and numbers of different 
types 

• added R glue code for various interfaces 
005 03/06/16 

• writeup and examples aplOuterProductO 

• writeup and examples aplRavelO 

• writeup and examples aplRankO 

• writeup and examples aplReduceO 

• writeup and examples aplReplicateO 

• writeup and examples aplReshapeO 

• added aplRcpp.hpp to code section 

006 03/06/16 

• writeup and examples aplGetO and aplSetO 
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• writeup and examples aplShapeO 


007 03/07/16 

• writeup and examples aplRotateO 

• fixed small bug, changed some variable names in aplRotateO 

• removed some unused variables from aplDotC.c 

• fixed bug in APLSCAN in aplDotCall.c 

008 03/07/16 

• writeup and examples aplScanO 

• correct small bug in aplReplicateO 

• writeup and examples aplSelectO 

• correct bug in aplTakeO 

• writeup and examples aplTakeO 

• writeup and examples aplTransposeO 

009 03/08/16 

• number_sections in pdf version 

• various typos corrected 
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